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Abstract 

We compute spin-dependent decay matrix elements for S-wave charmonium and bottomonium in 
lattice nonrelativistic quantum chromodynamics (NRQCD). Particular emphasis is placed upon the 
color-octet matrix elements, since the corresponding production matrix elements are expected to 
appear in the dominant contributions to the production cross sections at large transverse momenta. 
We use three slightly different versions of the heavy-quark lattice Green’s functions in order to 
minimize the contributions that scale as powers of the ultraviolet cutoff. The lattice matrix elements 
that we calculate obey the hierarchy that is suggested by the velocity-scaling rules of NRQCD. 
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I. INTRODUCTION 


The production rates of the and the T at large transverse momentum (p^) at the 
Fermilab Tevatron provide important tests of our understanding of heavy-quarkonium sys¬ 
tems and of quantum chromodynamics (QCD) itself. The dominant mechanism for both 
JI'll) and T production at large pr is expected to be the fragmentation of a gluon into a 
heavy-quark {QQ) pair JQ|. Because the heavy-quark mass m provides a large momentum 
scale, in addition to pr, the production of this heavy-quark pair can be calculated pertur- 
batively. The evolution of the QQ pair into the JM or the T is a nonperturbative process. 
In the color-singlet model P, 0, S, S fl Is, S 3 till; 3,13|, h is assumed that the 

QQ pair is in a color-singlet state. That assumption leads to predictions that are more than 

uarkonium 


an order of magnitude below the Tevatron data A less ad hoc approach to c 
production is based on the effective field theory nonrelativistic QCD (NRQCD) 


17 


mm- 


This approach is known as NRQCD factorization. The predictions of NRQCD factorization 
£t the Tevatron data and indicate that the dominant production mechanism at large pr 
proceeds through the fragmentation of a gluon into a QQ pair in a color-octet state. 

The effective field theory NRQCD separates long-distance, nonperturbative quarkonium 
dynamics (p ^ mv) from short-distance, perturbatively calculable processes {p^ m). Here 
m is the heavy-quark mass, p is the magnitude of the relative three-momentum of the Q 
and Q in the quarkonium rest frame, and v is the typical relative velocity of the Q and Q 
in the quarkonium rest frame. If A is the ultraviolet (UV) cutoff of NRQCD, then physics 
with p < A ~ m is reproduced in the effective field theory, while physics with p > A is 
integrated out, but affects the coefficients of local interactions in the effective theory. The 
Lagrangian of NRQCD can be expanded in powers of the velocity v. To any hnite order in 
V only a limited number of interactions appears in this Lagrangian. Hence, it is useful in 
calculations for systems in which -C 1. We note that in bottomonium ~ 0.1, while in 
charmonium ~ 0.3. 

In the context of NRQCD, scaling rules can be deduced for the leading behavior of opera- 
tor matrix elements in the limit in which the heavy-quark velocity v approaches zero 19.1^. 
When these “u-scaling rules” are applied to the production of the J/'ip and the T at large 
Pti they predict that, in the nonperturbative evolution of the QQ pair into the quarkonium, 
the non-spin-flip interactions dominate over the spin-flip interactions, with corrections of 
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order . Hence, the J/'ip is predicted to take on most of the transverse polarization of 


the gluon 


2l|. However, the CDF data for the polarization 22| show decreasing transverse 




polarization with increasing px and disagree with the NRQCD prediction |23| in the largest 
Pt bin. 

The existing calculations of polarization at the Tevatron neglect spin-flip processes in 
the NRQCD matrix elements, under the assumption that, because of the n-scaling rules, 
the spin-flip contributions are relatively suppressed. However, the n-scaling rules predict 
only the leading power of n, not its coefficient. It is usually assumed in making estimates 
that the coefficients are of order unity, but it could happen that the coefficients of the 
spin-flip terms are anomalously large or that the coefficients of the non-spin-flip terms are 


anomalously small. It has also been sugges 
to be modihed in the case of charmonium 
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tha 
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the n-scaling rules themselves may need 
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28|. 


One would like to test the applicability of estimates based on the n-scaling rules by 
direct calculation of the relevant NRQCD operator matrix elements. At present, lattice 
QCD is the only technique that is available for calculating color-octet matrix elements. 
Unfortunately, it is not known how to formulate that calculation of color-octet production 
matrix elements in Euclidean lattice held theory. However, one does know how to calculate 
the corresponding decay matrix elements in lattice NRQCD 0 y, 13 ]. Since the decay 
matrix elements are predicted to obey the same velocity-scaling rules as their production 
counterparts, a test of u-scaling estimates in the context of decay matrix elements might 
shed some light on the validity of n-scaling estimates for production matrix elements. In 
this paper, we present lattice NRQCD calculations of the decay matrix elements that are 
related by crossing to those production matrix elements that are expected to dominate S- 
wave quarkonium production at large pr- Preliminary results of these calculations were 
presented in Ref. j3l|. 


We use the methods developed by Lepage et al 


. 0 


to compute heavy-quark Green’s 


functions in lattice NRQCD. From these, we construct quarkonium propagators and use the 


>. y 


methods of Refs. mm I to extract the spin-dependent color-octet matrix elements that 
contribute to the decays of the J/'0, //c, and rjh (F-wave quarkonia) in lattice NRQCD. 
We perform the lattice measurements on a set of 400 12^ x 24 quenched conhgurations at 
f3 = Q/g"^ = 5.7. This f3 value was chosen because it corresponds to a momentum space 
cutoff A ~ mi,. In the case of charmonium one would like to have A ~ me- However, for 
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lattices that are significantly coarser than those with (3 = 5.7, it is not clear, even with 
highly improved actions, that one is close enough to the continuum limit to make contact 
with perturbation theory. We therefore use (3 = 5.7 lattices for charmonium and study some 
of the effects of those contributions that diverge as powers of A/rric by varying the algorithm 
for calculating lattice heavy-quark Green’s functions. 

Our bottomonium measurements suggest that the NRQCD n-scaling rules are useful in 
determining which matrix elements are most important, and that those contributions that 
are suppressed as powers of v show an even greater suppression than one would expect from 
the n-scaling factors alone. For charmonium the situation is complicated by our inability 
to use A ~ rric, but there are indications that the n-scaling rules are also a good guide 
here. Calculation of the perturbative coefficients that relate lattice and continuum matrix 
elements could help to clarify this situation. 

In Sec. m we introduce the NRQCD Lagrangian through relative order NRQCD 
factorization of quarkonium production and decay rates, and the lattice implementations of 
NRQCD that we use. We discuss the calculations that we performed and present our results 
in Sec. uni Sec.im contains discussions and our conclusions. 

II. NRQCD 

A. Continuum NRQCD 

As we indicated in the Introduction, NRQCD is an effective field theory with a UV 
momentum-space cutoff A ~ m. It is useful in describing bound states of heavy quarks. In 
the case of QQ bound states (quarkonium), the terms in the effective Lagrangian can be 
classified according to their leading power behavior in v. where v is the typical heavy-quark 
(or antiquark) velocity in the quarkonium rest frame |l SL l2flj | . The terms of leading order in 
in the NRQCD Lagrangian density are just the Schrodinger Lagrangian density: 



( 1 ) 


where = dt + igA^, D = d — igA, -0 is the Pauli spinor field that annihilates a heavy 
quark, and y is the Pauli spinor field that creates a heavy antiquark. 

In order to reproduce QCD completely, we would need an infinite number of interactions. 
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For example, at next-to-leading order in we have 


^>Cbilinear = ^ 

+ ^ -g^-g^- D)V^ + x^(D ■ pE - ■ D)x] 

-|- X g'Fj — 5fE x iD) • cr-^ -|- 'x^iiY) x gfE — ^fE x iD) ■ crx] 

+ ^ • (T)ij - x^(c/B ■ <t)x] . (2) 


In practice, we work to a given precision in v. For the calculations presented in this paper, 
the contributions given above suffice. 

It has been conjectured that, at large transverse momentum pt, the inclusive quarkonium 
production cross section can be written in a factorized form (Ref. jl^h 

( 3 ) 

n 

The “short-distance coefficients” Fn{A) are essentially the partonic cross sections to make a 
QQ pair with a given set of quantum numbers convolved with parton distributions. These 
partonic cross sections can be calculated as an expansion in The multiply vacuum 
matrix elements of four-fermion operators of the form 


o" = ,yK..V> 15^ l^^ + X){H + X\j ( 4 ) 

K contains Pauli matrices, color matrices, and the covariant derivatives Dt = dt + igA^, D = 
d — igA.} The operator matrix elements contain all of the long-distance, nonperturbative 
physics and are, essentially, the probabilities for a QQ pair with a given set of quantum 
numbers to evolve into a heavy-quarkonium state. The matrix elements are universal, i.e., 
process independent. NRQCD predicts the leading scaling behavior of the matrix elements 
with V (Ref. 13 )• ^ consequence of these n-scaling rules, the sum over operator matrix 

elements can be regarded as an expansion in powers of v, where ~ 0.3 for charmonium 
and 0.1 for bottomonium. 


^ A recent study of certain two-loop contributions to quarkonium production |32| has revealed that, if 
factorization is to hold, then the matrix elements must be modified from the form given in Eq. 0 by 
the inclusion of lightlike eikonal lines that run from each of the QQ bilinears to the far future. It is not 
known if this modification preserves the factorized form in higher orders. 
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A similar factorization formula applies to inclusive quarkonium decays |l^ : 

r(ffLH) = ^ 




nin 

n W 


( 5 ) 


except that the matrix elements are now between qnarkoninm states, rather than vacunm 
states, and the four-fermion operators have the form KnXX^i^'n'4’- While the coeffi¬ 

cients /n(A) can be calculated perturbatively, the quarkonium matrix elements {H\On{A.)\H) 
are nonperturbative and can be calculated directly in lattice NRQCD. An important fea¬ 
ture of NRQCD factorization is that both quarkonium decay and production occur through 
color-octet, as well as color-singlet, QQ states. While the production matrix elements are 
the crossed versions of quarkonium decay matrix elements, only the color-singlet production 
and decay matrix elements are simply related. 

The NRQCD operator matrix elements in Eqs. © and depend explicitly on the 
cutoff A. In the physical decay and production rates, this cutoff dependence in the matrix 
elements is canceled by a corresponding cutoff dependence in the short-distance coefficients. 
However, we note that the n-scaling rules for NRQCD matrix elements are derived under 
that assumption that the UV cutoff A is of order mv. The matrix elements of NRQCD, 
in common with those of other effective held theories, contain contributions that diverge, 
in the limit A —> cxo, both as logarithms and powers of A/m. When A is larger than mv, 
these divergent contributions potentially spoil the n-scaling rules, with the power divergences 
being especially important numerically. 

As an effective held theory, NRQCD is expected to be valid up to values of A close to 
m Q, 1^ . For bottomonium, A mj, is large enough that the physics at momenta greater 
than A can be treated perturbatively. While this choice should not invalidate the n-scaling 
rules for UV convergent contributions to matrix elements, it remains to be seen whether 
the n-scaling rules remain valid for matrix elements whose leading contributions diverge as 
A —>• cx). The measurements described in this paper test this conjecture, as well as the v- 
scaling rules for A ~ mv. In lattice evaluations of the NRQCD matrix elements, the value of 
the cutoh is determined, in part, by the lattice spacing. As we shall see, the ehective cutoh 
in the matrix elements is also affected by the specihc forms of the lattice action and Green’s 
functions that are employed.^ This latter property allows one to control the effective cutoff in 


^ In the continuum phenomenology of quarkonium production, the short-distance coefficients are usually 
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NRQCD matrix elements without affecting the interactions that determine the quarkonium 
masses and wave functions. Such an approach is especially useful in the case of charmonium, 
for which the cutoff A = rric is too small to include all of the nonperturbative bound-state 
physics. 

B. Lattice NRQCD 

The first part of this subsection summarizes the lattice formulation of NRQCD that has 
been given by Lepage et al. 0. 

In order to produce a lattice formulation of NRQCD, one must first formulate NRQCD 
in Euclidean space. This can be accomplished by performing the following substitutions in 
the NRQCD Lagrangian: 

—it, 
idt, 
iDt, 

-iE. (6) 

The gauge fields are incorporated into unitary matrices U^{x), which, as usual, are defined 
on the links of the lattice. Covariant derivatives are replaced on the lattice by covariant finite 
differences, which are defined by 

= + i)-Ul{x- i)'ilj{x - ?)], 

= ]^\Ui{x)ii{x + i) - 2ip{x) + Ul{x- i)ii{x - ?)], 

- 

i 

i 

calculated in dimensional regularization. In dimensional regularization, power infrared divergences in the 
short-distance coefficients are set to zero. This implies that, in the corresponding continuum NRQCD 
matrix elements, the UV power-divergent contributions are removed order-by-order in perturbation theory. 
Therefore, if the lattice matrix elements are to be close in value to the continuum matrix elements, then 
the effective A in the lattice calculations must be of order the heavy-quark mass or less. 


Ai'0(x) 

AfV(2:) 

A(2) 

A(4) 


t —>■ 
dt ^ 
Dt ^ 

</> ^ 
E ^ 
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Here we have adopted the standard lattice convention of working in a system of nnits in 
which the lattice spacing a has been set to nnity. 

Tadpole improvement is implemented by making the replacement 


U„{^) 


Mo 


Q 


( 8 ) 


For the tadpole-improved 


with Mo chosen as the fonrth root of the average plaqnette 
action, we expect the pertnrbation series for qnantities at the scale A to converge well 
Therefore, we replace each of the coefficients Cj in 5£biiinear with its lowest-order valne, 
namely, nnity. 

The order-M^ lattice Hamiltonian is now 


Ho = 


A(2) 

2m 


(9) 


The order-M"^ corrections to this Hamiltonian are 

( A ( 2))2 


6H = 


8m^ 

9 


+ 


W 

8m2 


A ■ E - E ■ A) 


8m2 

AH) 


(A X E-E X A) - — cr- B 

^ ^ 2m 


(A(2))^ 


24m 16nm2 


( 10 ) 


All of the terms except the last two are simple discretizations of those in the continnnm 
expression for (5£biiinear [Eq. (j21)]- The second-to-last term is the order-a^ correction to 
the discretization of the operator in Hq. The last term is the order-a^ correction to 
the approximation exp(—i^o) ~ [1 — Hq/{ 2n)Y^, which is used below in computing the 
evolution of a heavy-quark Green’s function over one lattice time step. Here, n is an integer 
to be specihed below. We note that the form of 5H in Eq. m is not tadpole-improved 
correctly by the replacement (0. That is because, in the higher-order derivatives A^^^ and 
there are canceling factors WU = 1 that are, incorrectly, replaced with 1 /mq. The 
corrections to this replacement amount to constant shifts of the Hamiltonian. Nevertheless, 
we retain the form of 6H in Eq. m, with the replacement (jHI), because, as noted by the 
NRQCD collaboration, it makes the contributions from 6H to the evolution of the heavy- 
quark Green’s functions small. The smallness of the 6H contributions can be discerned from 
the fact that the spin-averaged “masses” of the quarkonia are little affected by the inclusion 
of these order-M*^ corrections. 
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We make use of two different forms of the heavy-quark Green’s functions, which are 
equivalent through order One is the form that was used in the early spectroscopy papers 


of the NRQCD collaboration |30|, |3^. In this form, for a source S{x)Stfl, the retarded 
Green’s function Gr{x,t) at positive time is given recursively by 


Gr{x,0) = S{x)6tfi, 


Hn 




ul 


1-^1 G,(x,0), 


Gr{x,t+ 1 ) = ( 1 


H, 


n 

2,7' 


2n 


{1- 6H)Grix,t). 


( 11 ) 


We call this form of the Green’s function the “nrqcd scheme.” A second form of the Green’s 
function is dehned for positive time by 


Gr{x,t + 1 ) = ( 1 


Grix,0) = S{x)6tfi, 

This second form was used by the NRQGD collaboration in their heavy-light meson cal- 


Gr{x,t). (12) 


culations 


35|. For our calculation, we modify it with the rule that we replace any factor 


(1 — 6H/2) that appears as a result of consecutive time-evolution steps with {1 — 6H). This 
replacement results in an equivalent Green’s function to the order in v to which we work. 
We implement it because it ensures that the quarkonium spectra of the nrqcd scheme and 
this second scheme, which we call the “hybrid scheme” are identical. In both the nrqcd and 
hybrid schemes, we take Gr{x,t) = 0 for f < 0. In each scheme, we dehne a correspond¬ 
ing “advanced” Green’s function Ga, which satishes the same evolution equations as the 
hermitian conjugate of Gr, but with the boundary condition that it vanishes for f > 0. 

We note that the nrqcd and hybrid Green’s functions correspond to the same heavy-quark 
and antiquark propagation, except in the initial and hnal time slices. The hybrid scheme 
applies 6H to every time slice, including the initial and hnal time slices. In contrast, the 
nrqcd scheme does not apply 6H to the initial and hnal time slices, but applies 6H to all 
of the other time slices. In the operator matrix elements that we measure, the four-fermion 
operator is at the sink of the heavy-quark and antiquark propagators. Hence, in the nrqcd 
scheme, the interactions in 6H, which include those that change the spin, are turned oh for 
one time step on either side of the four-fermion operator. 

In order to see the ehect of the nrqcd scheme (and related schemes to be described later), 
let us initially ignore all of the interactions of the heavy-quark and antiquark with the gauge 
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field A, including those in Hq and in the time step on either side of the four-fermion 
operator. We call this variant the nrqcdx scheme. In the nrqcdx scheme, the vertices at which 
A interacts with the heavy-quark vanish on the time slice that contains the four-fermion 
vertex. Since this condition is a restriction on the temporal distance between the four- 
fermion vertex and the nearest heavy-quark-gluon vertex, its effect can be incorporated into 
the heavy-quark propagator that connects these vertices.^ One simply sets this propagator to 
zero at zero temporal separation. Now let us examine how this affects the momentum-space 
Feynman rule for the propagator. The momentum-space free-held heavy-quark propagator 
is given by the Fourier transform of the temporal propagator exp[—k^t/(2m)], which has 
support only for t > 0: 


G'o(fco,k) = 5^e^*'^°*exp 
t=o 



(13) 


Here k is the heavy-quark or antiquark momentum, the upper (lower) sign is for the quark 
(antiquark), and we have approximated the lattice Hamiltonian by the continuum expres¬ 
sion. In the nrqcdx scheme, the vanishing of the propagator at zero temporal separation is 
implemented by removing the t = 0 contribution from the above sum. That is, one replaces 
a free-held heavy-quark propagator adjacent to the four-fermion vertex with 


Gi(fc„,k) = ^ 


t=l 


„=F*fco4exp ( 


= exp ( =Fffco - ^ ) Go(fco,k). 


(14) 


On the right-hand side of Eq. da, the oscillatory factor exp(±ifco) provides additional 
damping of the integration, while the Gaussian factor exp[—k^/(2m)] cuts oh the inte¬ 
gration over the spatial components of the momentum at |k| ~ ^/rn. In the absence of these 
factors, the cutoh is of the order of the maximum value of the components of the lattice 
momentum, he., vr. In the nrqcdx scheme, the factors in Eq. m cut oh diagrammatic 
loops that contain interactions in 6H, but only for those loops that involve the four-fermion 
vertex. We note that these factors appear twice in each loop—once for each of the two 
heavy-quark propagators in the loop that attach to the four-fermion vertex. 


Note that additional vertices along the heavy-quark Green’s function cannot enter the time slice of the 


four-fermion interaction because the heavy-quark propagators have a definite time ordering. 
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As we have mentioned, the NRQCD matrix elements that we measure contain contribu¬ 
tions to that grow as powers of the lattice cutoff. These arise in the lowest nontrivial order 
in perturbation theory from loops involving the four-fermion vertex. Hence, the change in 
the effective cutoff for such loops that is provided by the nrqcdx scheme helps to control 
the numerical size of the power-behaved contributions in our simulations. We emphasize 
that the change in the effective cutoff affects only the interactions that renormalize the 
four-fermion operator. It has no effect, for example, on the interactions that produce the 
quarkonium masses or wave functions. 

It is useful to examine the effect of the nrqcdx scheme in coordinate space. If the four- 
fermion operator is at t = 0, then first interactions after the four-fermion interaction do 
not occur until t = 1, at which time the heavy-quark Green’s function is just the free-held 
Green’s function 

(15) 

Hence, in coordinate space, the effect of the nrqcdx scheme is to smear out the four-fermion 
interaction spatially over a distance of order 1/y/m. 

Now let us return to the effects of including, in the time slices that are adjacent to the 
four-fermion interaction, the interactions with the gauge held A that are contained in Hq and 
Ui- This discussion is relevant to the nrqcd scheme and to the coulomb scheme, which we 
introduce later. By ignoring the ehects of the Goulomb gluon held Aq = 0, we have ignored 
the ehects of the potential in the Schrodinger equation. This is a reasonable approximation 
until the separation of the QQ pair is of the order of the typical QQ separation in the bound 
state, at which point one can take the spatial smearing to be given roughly by the size of 
the QQ bound state. Since we work in the Goulomb gauge, the ehects of the heavy-quark 
interactions with the spatial gauge held A are subleading in v relative to the ehects of the 
potential. 
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III. CALCULATIONS AND RESULTS 


We are interested in matrix elements between quarkonium states of operators of the form 


and 


Oi (^ 5 * 0 , 0 ) = (16a) 

Oi f = 't/;^a^xx^^±'(p, (16b) 

Oi f 5'i,o) = V’'^cr3XX^o-3^, (16c) 

Cls C^o.o) = (17a) 

Cls (17b) 

Os (^5*1,0) = ^/’VsT^xx^o'sT'^, ( 17 c) 


where the subscripts 1 and 8 indicate color-singlet and color-octet operators, respectively, 
and S denotes an S'-wave operator. In s is the total spin quantum number, m is the 

quantum number of the component of the spin along the quantization axis, T“ is an SU(3) 
color matrix in the fundamental representation satisfying Tr(T“T^) = (1/2)5“^, and the ex’s 
are Pauli spin matrices. In Eqs. (El), there is an implied sum over the color index a. The 
matrix elements of the operators in Eqs. (HED and ED in quarkonium states are those that 
appear in the NRQCD factorization formula for decays dSD- We note that, according to the 
n-scaling rules of NRQCD, the matrix elements of the color-octet operators are suppressed 
by at least relative to the leading color-singlet matrix elements. 

We are interested in the matrix elements of the operators in Eqs. (IWH) between both 
the spin-singlet (pseudoscalar) and spin-triplet (vector) S'-wave quarkonium states. What 
we measure is the ratio of octet to singlet matrix elements 




/2si+l q 

\ ^Si.rrii 

1^8 

) 

|2si+lc< \ 

1 ^Si,mi/ 

/2si+lC 1 


mi/vs 


(18) 


where the initial (hnal) spins and 2 ;-components of spin are Si (s/) and m* (m/), respectively. 
The subscript VS indicates that we have used the vacuum-saturation approximation [l^ for 
the denominator. That is, we replace the sum over intermediate states in the operator with 
the vacuum state. We do, in fact, measure the corresponding ratios for the case of color- 
singlet operators in the numerator, and we hnd that the spin-diagonal ratios are close to 
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one, indicating that onr use of the vacuum-saturation approximation is valid and that the 
off-diagonal matrix elements are all small, as is predicted by v scaling. 

Now let us describe how we measure R on the lattice. First, we create a QQ pair, using 
a source on a time slice t. We propagate this pair to a second time slice t' > t, using the 
equations for the retarded Green’s functions given in Sec. IHI Then we annihilate the QQ pair 
at a point, in the spin and color state that corresponds to the operator of interest. We then 
re-create the QQ pair at the same point and propagate this pair to time slice t" > t', where it 
is annihilated at a sink. In practice we use sinks on the time slice t" as sources and propagate 
the QQ pair back to t' by using the “advanced” Green’s functions, conjugating these to give 
us the Green’s functions that we need. We call the quantity that we have just described 
the “lattice-matrix-element precursor.” For T = t' — t and T' = t" — t' sufficiently large, we 
annihilate the QQ pair from an almost pure quarkonium ground-state wave function. The 
quarkonium propagator falls as exp(—i?T), where E is the energy of this ground state. By 
computing a ratio of color-octet to color-singlet lattice-matrix-element precursors, we cancel 
this exponential falloff, as well as the amplitude factors that are associated with the overlaps 
between the sources and sinks and the quarkonium wave function. The result, in the limit 
of large T and T', is precisely the lattice version of the ratio R. 

In our calculations, we use one of two stochastic sources on the initial time slice to 
generate our retarded Green’s functions. The hrst of these sources consists of a complex 
random number that is uniformly distributed in U{1) at each site. The second of these 
sources consists of a Gaussian smearing of this U{1) random source. By choosing the Q and 
Q Green’s functions to have conjugate U{1) sources, we obtain a stochastic estimator of a 
point source for the QQ pair at each site on the initial time slice. By choosing the the U{1) 
source for the Q and the conjugate of the Gaussian-smeared U{1) source for the Q (or vice 
versa), we obtain a stochastic estimator of a Gaussian source for the QQ pair at each site 
on the initial time slice. The advanced Green’s functions are treated similarly. 

We calculate the quarkonium propagators and matrix elements on each of 400 equilibrated 
12^ X 24 quenched lattices at 0 = 5.7. We use the parameters that were determined by the 


NRQGD collaboration 


3G 


M 


36|. For b quarks, we use m = 3.15, and for c quarks we use 


m = 0.8. We take Uq = 0.860846184, which we obtained from our own measurements of 
the average plaquette. We choose for the width of our Gaussian source 2.5 lattice units, for 
both charmonium and bottomonium. Fig. ^ shows effective quarkonium “mass” (energy) 
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plots for point and Gaussian-smeared sources as functions of the time separation T between 
the source and the sink. These plots indicate that the Gaussian width that we take is 
reasonable, as it leads to an early approach to the asymptotic value of the energy. The hts 
to the effective energy in the nrqcd scheme with the extended source/sink give an energy 
of Eq = 0.6165(6), with a conhdence level of 78% for a £t over the T(T')-range 8-19, and a 
J/'ip energy of Ei = 0.6938(7), with a conhdence level of 80% for a fit over the T(T')-range 
10-22. The corresponding hts for bottomonium give an r]b energy of Eq = 0.5065(3), with 
a conhdence level of 60% over the T(T')-range 8-22, and an T energy of Ei = 0.5215(3), 
with a conhdence level of 59% for a ht over the T(T')-range 5-22. These results are to be 
compared with those of the NRQGD collaboration, which obtained Eq = 0.6182(7) for the 
r]c, El = 0.697(1) for the J/ip, Eq = 0.5029(5) for the rjb, and Ei = 0.5186(6) for the T. 

Now let us turn to the matrix elements. Here we recall that the decay matrix elements 
that we measure are related to the production matrix elements that are of primary interest 
through the crossing of the quarkonium from the intermediate state to the initial state. In 
the production of the J/p) oi the T by gluon fragmentation—the dominant process at large 
Pt —the gluon fragments into a QQ pair in a triplet spin state with transverse polarization. 
Hence, we are most interested in the decay matrix elements in which the annihilated QQ 
pair is transversely polarized. 

We measure ratios R of color-octet matrix elements to color-singlet matrix elements 
dehned in Eq. Ha. For these measurements we use both the hybrid and the nrqcd schemes 
for calculating the required heavy-quark Green’s functions. 

Let us hrst examine the bottomonium matrix elements. Here the spin-singlet quarkonium 
is the rjb and the spin-triplet quarkonium is the T. We ht the measured ratios over a range 
of T and T' values with T -I- T' < 24 in order to maximize the information that we extract 
from sometimes noisy “data.” Our measurements are shown in Table HI 

The errors presented in Table H] are statistical only. For the matrix elements that connect 
the singlet and triplet spin states and, in the hybrid scheme, for the up —> down matrix 
element, the signal is excellent, and the systematic errors associated with our choices of hts 
are probably less than the statistical errors. Such estimates of htting errors are obtained by 
examining the hts over a number of choices of htting ranges in T and T' (132 choices for 
each matrix element) and taking into account the range of htted values and statistical errors 
for the matrix element of interest. Where possible, we consider only those hts for which 
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TABLE I: Ratios R of bottomonium color-octet matrix elements to color-singlet matrix elements, 
as defined in Eq. m- The label singlet denotes a state 5^o,o? triplet denotes a state 
“up” denotes a state “down” denotes a state and “longitudinal” (long.) denotes a 

state The spin transition i —> / is from an initial state i to an annihilated state / in the 

numerator of R. In the case of triplet states, we average over rrii in initial states and sum over m/ 
in annihilated states. The labels “nrqcd” and “hybrid” refer to the methods for calculating the 
heavy-quark lattice NRQCD Green’s functions. The column labeled “u scaling” gives the u-scaling 
factors of the various ratios R, along with their numerical values. In evaluating the u-scaling 
factors, we take = 0.1 for bottomonium. 


spin transition 

hybrid 

nrqcd 

V i 

scaling 

singlet - 

-> triplet 

7.169(6) X 10-3 

2.72(4) X 10-^ 

u3/(2iVe) 

« 5.3 X 10-3 

triplet - 

singlet 

2.414(3) X 10-3 

9.0(1) X 10-3 

vV{2Nc) 

ss 5.3 X 10-3 

singlet - 

-> singlet 

6.1(5) X 10-3 

6.5(5) X 10-3 

uV(2W) 

« 1.7 X 10-3 

triplet - 

triplet 

8.1(6) X 10-3 

6.9(5) X 10-3 

vV{2N,) 

ss 1.7 X 10-3 

up - 

up 

7.3(6) X 10-3 

6.9(6) X 10-3 

vVm,) 

ss 1.7 X 10-3 

long. —> 1 

transverse 

2.7(3) X 10-3 

1 - 2 X 10-3 

vV{2N,) 

« 1.7 X 10-^ 

down 

^ up 

5.53(2) X 10-3 

< 5 X 10-3 

vV{2N,) 

« 1.7 X 10-^ 


T and T' are large enough that the value of the matrix element appears to have reached a 
stable plateau. In the remainder of this paper, when we refer to systematic errors, we mean 
systematic errors associated with our choices of fits, unless we explicitly state otherwise. 

For the three diagonal matrix elements in both schemes, in which the signal is noisy 
but still substantial, the systematic errors could be as large as 20%. The signal for the 
longitudinal —> transverse matrix element in the hybrid scheme is sufficiently weak that the 
systematic errors could be as much as 50%. Finally, the signals for the two triplet spin-flip 
matrix elements in the nrqcd scheme are so noisy that we are only willing to place bounds 
on their values. 

Now let us compare our results with expectations from the n-scaling rules of NRQCD. 
The n-scaling factors of the ratios R are given in the last column of Table d For purposes 


of using the n-scaling factors to estimate the sizes of matrix elements, we include 
factor l/{2Nc) that arises in the free QQ matrix elements, as suggested in Ref. 


he color 
37|. We 
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FIG. 1; The effective quarkonium “mass” (energy) as measured with point and Gaussian-smeared 
(extended) sources for the J/'0 and the T. E is the effective energy, and T is the time separation 
between the source and the sink. The effective energies shown are measured in the nrqcd updating 
scheme that is described in Sec. HH The plots in the hybrid scheme and for the rjc and the r]i, are 
similar. 
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evaluate the u-scaling factors by taking = 0.1 for bottomonium. Such u-scaling estimates 
are based on the assumption that the coefficients in the expansions of the matrix elements 
in powers of v are of order unity. It is that assumption that we wish to test by our explicit 
calculation. Of course, if it turns out that the coefficients are signihcantly greater than order 
unity, then the v expansion will be of little use. 

As we have mentioned, the NRQCD matrix elements that we measure contain contribu¬ 
tions that diverge as powers of A in the limit A —> cx). These power-divergent contributions 
potentially violate the n-scaling rules, which were derived under the assumption that A is 
of order mv. For the particular matrix elements that we measure, a perturbative analy¬ 
sis shows that the leading contributions to the singlet—triplet matrix elements diverge as 
[Q;5(A)/7r](A/mfe)^, the diagonal contributions diverge as [Q;s(A)/7r]^(A/mf,)^ log^(A/mfe), and 
the triplet spin-flip contributions diverge as [Q;s(A)/7r]^(A/mfe)^. Unless A is not much larger 
than mv, these power-behaved contributions may lead to signihcant numerical violations of 
the n-scaling rules. As we mentioned earlier, we wish to test not only whether the n-scaling 
rules hold for cutoffs A ~ mbV, but also whether they continue to hold for cutoffs A ~ mb- 
Therefore, in testing the n-scaling rules in the bottomonium system, we choose A/mb ~ 1- 
As we stated earlier, at /? = 5.7, the lattice momentum cutoff itself is close to the input 
6-quark mass. Hence, in the case of bottomonium, no special choice of heavy-quark Green’s 
functions is required in order to control power-behaved contributions, and we focus on the 
hybrid scheme. 

In the hybrid scheme, we note that the singlet —triplet matrix elements are the largest 
color-octet matrix elements, as is predicted by the u-scaling estimates. The u-scaling esti¬ 
mates also predict their magnitudes correctly. The heavy-quark spin symmetry of NRQCD 
predicts that the singlet —triplet matrix element should be a factor 3 larger than the 
triplet —singlet matrix element, with corrections to that relation of order This pre¬ 
diction is borne out by our measurements. The diagonal matrix elements are suppressed 
relative to the singlet —triplet and triplet —singlet matrix elements, as expected, but 
by considerably more than would be expected from the u-scaling factors alone, suggesting 
that their coefficients in the v expansion are small. The triplet spin-flip matrix elements are 
smaller still, as is expected from the u-scaling estimates, but, again, they are much smaller 
than would be expected from the u-scaling factors alone. The suppression of the triplet 
spin-flip matrix elements relative to the diagonal matrix elements is approximately what 
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would be expected from the u-scaling factors. 

The effects of a change in the effective cutoff A can be seen when we compare the hybrid- 
scheme and nrqcd-scheme values for R in Table HJ The nrqcd scheme results in a smaller 
effective cutoff for the interactions in 5H, which include all of the spin-dependent interac¬ 
tions. We see that the nrqcd-scheme singlet —triplet and triplet —>■ singlet matrix elements 
are suppressed by a factor of about 26 relative to the corresponding hybrid-scheme matrix 
elements. The triplet spin-flip matrix elements also appear to be suppressed. The diagonal 
matrix elements are virtually unchanged. That is expected because the leading contribution 
to the diagonal matrix elements comes from whose effective cutoff is not changed in 
going from the hybrid scheme to the nrqcd scheme. 

Now let us turn to the charmonium case. For charmonium, the heavy-quark mass in 
lattice units is only 0.8, while, as usual, the components of the lattice momentum can range 
up to TT in magnitude. Hence, we expect that the lowering of the effective cutoff A that 
is provided by using the nrqcd scheme will be helpful in reducing the effects of the power- 
behaved contributions. In the case of charmonium, we will also consider an additional 
scheme for computing the heavy-quark Green’s functions in which the links Ui{x) are set 
to unity on the time slice that is associated with the four-fermion operator (the time slice 
in which the QQ pair is annihilated and re-created). In other respects, this scheme, which 
we call the “coulomb scheme” is identical to the nrqcd scheme. In the coulomb scheme, we 
are neglecting the interactions of the heavy quark with the helds A, on the QQ'^iiiiihilation 
time slice. Since we work in the Coulomb gauge, these interactions are subleading in The 
effect of the coulomb scheme, relative to the nrqcd scheme, is to lower the effective cutoff 
on Hq so that it is the same as the effective cutoff on 5H . Hence, we expect the diagonal 
matrix elements to show the effects of a reduced cutoff. 

The values of R in the hybrid, nrqcd, and coulomb schemes, along with the n-scaling 
factors, are given in Table El Again, the quoted errors are statistical only. In the hybrid 
scheme, the signal for the triplet —> singlet, singlet —> triplet, and up —down matrix 
elements is robust, and we feel conhdent that the systematic errors are less than the statistical 
errors. The signal for the singlet—*• singlet matrix element is noisier, but still quite good, so 
that we feel that the systematic error is probably no worse than the statistical error. The 
signals for the up —> up and transverse —longitudinal matrix elements are even noisier, 
but still substantial. In these cases, the systematic errors might be as large as 10%. In the 
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TABLE II: Ratios R of charmonium color-octet matrix elements to color-singlet matrix elements, 
as defined in Eq. m- The labels are as in Tabled except that “coulomb” refers to an additional 
method for calculating the heavy-quark lattice NRQCD Green’s functions, as described in the text. 
Eor purposes of numerical estimates of the u-scaling factors, we take = 0.3 for charmonium. 


spin transition 

hybrid 

nrqcd 

coulomb 

V f 

scaling 

singlet - 

triplet 

6.397(8) X 10-2 

2.90(3) X 10-3 

2.84(4) X 10-3 

vVm,) 

Ri 2.7 X 10-2 

triplet - 

-> singlet 

2.938(7) X 10-2 

1.13(2) X 10-3 

1.06(1) X 10-3 

vV{2N,) 

Rs 2.7 X 10-2 

singlet - 

-> singlet 

5.03(9) X 10-^ 

9.7(2) X 10-^ 

3.6(3) X 10-^ 

vV{2N,) 

Rs 1.5 X 10-2 

triplet - 

-> triplet 

1.57(2) X 10-3 

1.016(8) X 10-3 

4.0(4) X 10-^ 

vV{2N,) 

Rs 1.5 X 10-2 

up - 

up 

4.57(6) X 10-^ 

1.019(8) X 10-3 

3.9(4) X 10-^ 

vV{2N,) 

RS 1.5 X 10-2 

long. —> 1 

transverse 

2.82(4) X 10-^ 

2.8(7) X 10-6 

— 

vV{2Nc) 

Ri 4.5 X 10-3 

down 

^ up 

8.48(5) X 10-^ 

1.4(2) X 10-6 

1.4(3) X 10-6 

vVi2Nc) 

RS 4.5 X 10-3 


nrqcd scheme, the signals for the triplet —>• singlet, singlet —triplet, and singlet —singlet 
matrix elements are as clean as those in the hybrid scheme. The signals for the two diagonal 
triplet matrix elements show some rise as T and T' are increased, and the systematic error 
could be as large as 10%. The signals for the two triplet spin-flip matrix elements are so 
noisy that the values given should be considered only to be order-of-magnitude estimates. 
In the coulomb scheme, the signals for all of the matrix elements are noisier than in either 
of the other two schemes. The triplet ^ singlet and singlet —> triplet elements probably 
have systematic errors that are comparable with their statistical errors. The diagonal matrix 
elements could have systematic errors as large as 20%, while the signals for the triplet spin- 
flip matrix elements are so noisy that we can only say that these matrix elements are small. 
In Figure 01 we show examples of the signals for the charmonium ratios R as functions of 
T = T'. The upper plot shows R for the singlet —> triplet transition in the hybrid scheme, 
which has one of the clearest signals that we have seen. The lower plot shows R for the 
longitudinal —>■ transverse transition in the coulomb scheme, which has one of the noisiest 
signals that we have seen. 

Now let us compare these results with the predictions from the n-scaling factors, which are 
shown in the last column of Tableim As in the bottomonium case, we include the color factor 
l/{2Nc) that arises in the free QQ matrix elements (Ref. 0 ), The 
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FIG. 2: Values of ratios R for charmonium as a function of T = T'. The upper plot shows R 
for the singlet —> triplet spin transition in the hybrid scheme. The lower plot shows R for the 
longitudinal —> transverse spin transition in the coulomb scheme. 
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factors are the same as those for bottomonium, except that we evaluate them by setting 
= 0.3. Again, as in the case of bottomonium, the triplet —> singlet and singlet —triplet 
matrix elements in the hybrid scheme show good agreement with the n-scaling predictions. 
The other matrix elements are somewhat smaller than would be expected from the n-scaling 
predictions, perhaps indicating that their coefficients in the v expansion are small. However, 
here we notice that the triplet spin-flip and non-spin-flip matrix elements are of similar size. 
As we have mentioned, we know that the effective cutoff in the charmonium case is larger 
than m and that the spin-flip matrix elements diverge as [Q;5(A)/7r]^(A/mc)"^ at large A, while 
the non-spin-ffip matrix elements diverge only as [Q;s(A)/7r]^(A/mc)^ log^(A/mc). Therefore, 
it is not surprising that the triplet spin-flip matrix elements appear to be anomalously large 
in comparison with the non-spin-flip matrix elements. In the nrqcd scheme, which effectively 
lowers the cutoff on the interactions in 5H, we see a considerable reduction in the size of 
the triplet spin-flip matrix elements, as well as in the triplet —singlet and singlet —triplet 
matrix elements, which diverge as [Q;s(A)/7r](A/mc)^. However, the sizes of the diagonal 
matrix elements are not reduced, but merely made closer in value, as would be expected 
when all off-diagonal matrix elements are small. This indicates that these diagonal matrix 
elements get sizable contributions from the interactions in Hq, whose effective cutoff is not 
reduced in the nrqcd scheme. In the coulomb scheme, in which the effective cutoff of the 
interactions in Hq is also reduced, we see that the diagonal matrix elements are reduced 
in size relative to those in the nrqcd scheme. The triplet —singlet and singlet —triplet 
matrix elements are almost identical in size to their nrqcd-scheme counterparts because 
the effective cutoff on the interactions in 5H is the same (to leading order in v) in both 
schemes. In addition, to the extent that we can infer anything from the very noisy signals 
for the triplet spin-flip matrix elements in the coulomb scheme, the values of these matrix 
elements are consistent with their remaining unchanged in going from the nrqcd scheme to 
the coulomb scheme. We can therefore conclude that, were we able to carry out the lattice 
computation directly at A ~ me, the size ordering of matrix elements would be as predicted 
by the n-scaling estimates. However, the suppression of the matrix elements of higher order 
in V would be even more than that which is predicted by those estimates, indicating that 
the those matrix elements have relatively small coefficients in the v expansion. We conclude 
that the v expansion is a useful tool for establishing a hierarchy of charmonium decay matrix 
elements. 
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IV. SUMMARY AND CONCLUSIONS 


We have calculated, in lattice NRQCD, color-octet, spin-dependent matrix elements that 
appear in the NRQCD factorization expression © for decays of bottomonium and char- 
monium. The lattice action we used is correct through order and contains the standard 
tadpole improvement, as well as order-a^ improvements to the terms of leading order in 
V. The decay matrix elements that we calculated are related by crossing symmetry to the 
production matrix elements that appear in the dominant contributions at large pt to the 
production of the T and the J/'0- In our calculations, we made use of various forms of 
the lattice heavy-quark Green’s functions, which are equivalent to the order in v to which 
we work, but which allow us to vary the effective momentum cutoff for interactions of the 
heavy-quark helds with the gauge helds. Our goals were to test the accuracy of estimates 
of the sizes of matrix elements that are based on the u-scaling rules of NRQCD and to test 
the convergence of the v expansion. 

In the case of bottomonium, we found that estimates that are based on the u-scaling 
factors are reasonable for the singlet —> triplet and triplet —> singlet matrix elements. The 
diagonal singlet —singlet and triplet —triplet matrix elements are suppressed relative to 
these leading-order matrix elements, but by even more than one would expect from the 
u-scaling factors alone. The suppression of the triplet spin-flip matrix elements relative to 
the triplet non-spin-flip matrix elements is consistent with predictions that are based on the 
u-scaling factors. In all cases, the u-scaling predictions correctly indicate the hierarchy of 
matrix elements. 

In the case of charmonium, the situation is made less clear by the fact that, in order to 
describe the quarkonium physics on the lattice, we are forced to work with lattice spacings 
that correspond to UV momentum cutoffs that are considerably larger than the heavy-quark 
mass. Consequently, we expect the matrix-element calculations to contain large contribu¬ 
tions that grow as a power of the cutoff and that violate the u-scaling rules of NRQCD. By 
using various forms of the heavy-quark Green’s functions, we are able to vary the effective 
UV cutoffs of the matrix elements. The indications from our studies of the effects of varying 
these cutoffs are that the u-scaling rules are a good guide as to which matrix elements are im¬ 
portant, but that the matrix elements are smaller than one would expect from the u-scaling 
factors alone. Again, the u-scaling predictions correctly indicate the hierarchy of matrix 
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elements. Of course, the cutoff dependences of the matrix elements are canceled in physical 
quantities, such as the quarkonium production and decay rates, by corresponding cutoff de¬ 
pendences in the accompanying short-distance coefficients. A more complete understanding 
of the effects of varying the UV cutoffs could be obtained by studying the cutoff dependences 
of the matrix elements in perturbation theory. For example, one could calculate the match¬ 
ing coefficients between lattice matrix elements and continuum (dimensionally-regulated) 
matrix elements. The matching coefficients are infrared safe and, hence, are amenable to a 
perturbative treatment. Such a calculation involves both one- and two-loop contributions 
in lattice and continuum perturbation theory, and is beyond the scope of this paper. 

Phenomenological color-octet production matrix elements for the triplet —triplet tran¬ 
sition, normalized against the color-singlet matrix elements, as in R in Tables H] and IHl he 
in the ranges 5.1-16 x 10“^ for the T(IS') and 1.9-24 x 10“^ for the J/'ijj (Ref. j^). The 
lattice decay matrix elements for the triplet —triplet transition are somewhat smaller than 
the phenomenological range of production matrix elements for the J/V’ and considerably 
smaller than the phenomenological range of production matrix elements for the T. How¬ 
ever, effects from multiple-gluon radiation, which may decrease the phenomenological values 
of the color-octet matrix elements, have not yet been included in the analyses of the T ma¬ 
trix elements. Furthermore, decay matrix elements need not be equal to the corresponding 
production matrix elements, even though they have the same u-scaling behavior. We also 
note that lattice matrix elements differ from continuum matrix elements. We do not expect 
the differences to be large for tadpole-improved lattice calculations, provided that the lattice 
cutoff is of order or smaller than the heavy-quark mass. However, that situation does not 
hold in our charmonium calculations. 

The charmonium color-octet singlet —> triplet decay matrix element is relatively large, 
which suggests that the color-octet triplet —singlet production matrix element may also 
be large. This would imply that, at large pt at the Tevatron, where the color-octet, spin- 
triplet process dominates F-wave quarkonium production, the rjc production rate may be 
comparable to the J/ijj production rate. Once the size of the effective UV cutoff is reduced, 
the longitudinal —> transverse decay matrix element is small relative to the up —> up decay 
matrix element. This suggests that a similar hierarchy may hold for the transverse —>■ lon¬ 
gitudinal and up ^ up production matrix elements, which would support the prediction of 
large transverse polarization at large px at the Tevatron. 
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